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Abstract 


The interaction of light with matter has triggered the interest of scientists for a long time. 
The area of plasmonics emerges in this context through the interaction of light with 
valence electrons in metals. The random phase approximation in the long wavelength 
limit is used for analytical investigation of plasmons in three-dimensional metals, in a 
two-dimensional electron gas, and finally in the most famous two-dimensional semi- 
metal, namely graphene. We show that plasmons in bulk metals as well as in a two- 
dimensional electron gas originate from classical laws, whereas quantum effects appear 
as non-local corrections. On the other hand, graphene plasmons are purely quantum 
modes, and thus, they would not exist in a “classical world.” Furthermore, under certain 
circumstances, light is able to couple with plasmons on metallic surfaces, forming a 
surface plasmon polariton, which is very important in nanoplasmonics due to its 
subwavelength nature. In addition, we outline two applications that complete our 
theoretical investigation. First, we examine how the presence of gain (active) dielectrics 
affects surface plasmon polariton properties and we find that there is a gain value for 
which the metallic losses are completely eliminated resulting in lossless plasmon prop- 
agation. Second, we combine monolayers of graphene in a periodic order and construct 
a plasmonic metamaterial that provides tunable wave propagation properties, such as 
epsilon-near-zero behavior, normal, and negative refraction. 


Keywords: random phase approximation, graphene, gain dielectrics, plasmonic 
metamaterial 


1. Introduction 


The interaction of light with matter has triggered the interest of scientists for a long time. The 
area of plasmonics emerges in this context through the interaction of light with electrons in 
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metals, while a plasmon is the quantum of the induced electronic collective oscillation. In 
three-dimensional (3D) metals as well as in a two-dimensional electron gas (2DEG), the 
plasmon arises classically through a depolarized electromagnetic field generated through 
Coulomb long-range interaction of valence electrons and crystal ions [1]. Under certain cir- 
cumstances, light is able to couple with plasmons on metallic surfaces, forming a surface 
plasmon polariton (SPP) [2-4]. The SPPs are very important in nanoplasmonics and 
nanodevices, due to their subwavelength nature, that is, because their spatial scale is smaller 
than that of corresponding free electromagnetic modes. In addition to classical plasmons, 
purely quantum plasmon modes exist in graphene, the famous two-dimensional (2D) semi- 
metal. Since we need the Dirac equation to describe the electronic structure of graphene, the 
resulting plasmons are purely quantum objects [5-8]. As a consequence, graphene is quite 
special from this point of view, possessing exceptional optical properties, such as ultra- 
subwavelength plasmons stemming from the specifics of the light-matter interaction [7-10]. 


In this chapter, we present basic properties of plasmons, both from a classical standpoint but 
also quantum mechanically using the random phase approximation approach. Plasmons in 3D 
metals as well as in 2DEG originate from classical laws, whereas quantum effects appear as 
non-local corrections [11-13]. In addition, we point out the fundamental differences between 
volume (bulk), surface, and two-dimensional plasmons. We show that graphene plasmons are 
a purely quantum phenomenon and that they would not exist in a “classical world.” We then 
outline two applications that complete our theoretical investigation. First, we examine how the 
presence of gain (active) dielectrics affects SPP properties and we find that there is a gain value 
for which the metallic losses are completely eliminated resulting in lossless SPP propagation 
[3]. Second, we combine monolayers of graphene in a periodic order and construct a plasmonic 
metamaterial that provides tunable wave propagation properties, such as epsilon-near-zero 
behavior, normal, and negative refraction [9]. 


2. Volume and surface plasmons in three-dimensional metals 


2.1. Free collective oscillations: plasmons 


Plasma is a medium with equal concentration of positive and negative charges, of which at 
least one charge type is mobile [1]. In a classical approach, metals are considered to form 
plasma made of ions and electrons. The latter are only the valence electrons that do not interact 
with each other forming an ideal negatively charged free electron gas [1, 14]. The positive ions, 
that is, atomic nuclei, are uniformly distributed forming a constant background of positive 
charge. The background positive charge is considered to be fixed in space, and as a result, it 
does not respond to any electronic fluctuation or any external field while the electron gas is 
free to move. In equilibrium, the electron density (plasma sea) is also distributed uniformly at 
any point preserving the overall electrical neutrality of the system. Metals support free and 
collective longitudinal charge oscillation with well-defined natural frequency, called the 
plasma frequency wp. The quanta of these charge oscillations are plasmons, that is, quasi- 
particles with energy Ep = hw ,, where h is the reduced Plank constant. 
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We assume a plasma model with electron (and ion) density n. A uniform charge imbalance ôn 
is established in the plasma by displacing uniformly a zone of electrons (e.g., a small slab in 
Cartesian coordinates) by a small distance x (Figure 1). The uniform displacement implies that 
all electrons oscillate in phase [2]; this is compatible with a long wavelength approximation 
(A,/a — œ, where A, is the plasmon wavelength and a is the crystal lattice constant); in this 
case, the associated wavenumber |q| (Figure 1(b)) is very small compared with Fermi wave- 
number kr, viz. q/kp — 0 [7]. Longitudinal oscillations including finite wave vector q will be 
taken into account later in the context of quantum mechanics. The immobilized ions form a 
constant charge density indicated by en, where e is the elementary charge. Let x(t) denote the 
position of the displaced electronic slab at time t with charge density given by —edn(t). Due to 
the electron displacement, an excess positive charge density is created that is equal to edn(t), 
which in equilibrium, 6n = 0, reduces to zero. Accordingly, an electric field is generated and 
interacts with the positive background via Coulomb interaction, forcing the electron cloud to 
move as a whole with respect to the immobilized ions, forming an electron density oscillation, 
that is, the plasma oscillation. The polarized electric field is determined by the first Maxwell 
equation as 


V -E=4ne6n, (1) 


in CGS units.’ The displacement x(t) in the electronic gas produces an electric current density 


J = -e(n + 6n)x = — enx (since 6n/n — 0), related to the electron charge density via the conti- 
nuity equation V - J = —ed,6n. After integration in time, we obtain 
én =nV-x (2) 


Combining Eqs. (1) and (2), we find the electric field that is induced by the electron charge 
displacement, that is, 


(a) 


Figure 1. (a) A charge displacement is established by displacing uniformly a slab of electrons at a small distance x, 
creating a polarized electric field in the solid. (b) A plasma longitudinal oscillation electric field in the bulk of a solid. 
The arrow indicates the direction of displacement of electrons and of the wavevector q, while the double-faced arrow 
shows the plasmon wavelength Ap. 


‘For SI units, we make the substitution 1/¢9 = 47. 
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E = 47enx. (3) 


Newtonian mechanics states that an electron with mass m in an electric field E obeys the 
equation mx = —eE, yielding finally the equation of motion 


mx + 4re*nx = 0, (4) 


indicating that electrons form a collective oscillation with plasma frequency 


w,(0) = 4/ pe (5) 


where w,(0) = w,(q = 0). The energy Ep = hw, is the minimum energy necessary for exciting a 


plasmon. Typical values of plasmon energy E, at metallic densities are in the range of 2 — 20 eV. 


Having shown that an electron gas supports free and collective oscillation modes, we proceed 
to investigate the dynamical dielectric function ¢(q,w) of the free electron gas. The dielectric 
function is the response of the electronic gas to an external electric field and determines the 
electronic properties of the solid [1, 11, 15]. We consider an electrically neutral homogeneous 
electronic gas and introduce a weak space-time-varying external charge density p,,,(x,t) [14]. 
Our goal is to investigate the longitudinal response of the system as a result of the external 
perturbation. In free space, the external charge density produces an electric displacement field 
D(x,t) determined by the divergence relation V -D = 47p,,,. Moreover, the system responds 
and generates additional charges (induced charges) with density p,,(x,t) creating a polariza- 
tion field P(x,t) defined by the expression V - P = —p,,,4 [1]. Because of the polarization, the 
total charge density inside the electron gas will be Pot = Pext + Ping; leading to the screened 
electric field E, determined by V - E = 47p,,,;. The fundamental relation D = E + 47P is derived 
after combining the aforementioned field equations. 


The dielectric function is introduced as the linear optical response of the system. According to 
the linear response theory and taking into account the non-locality in time and space [2, 14], 
the total field depends linearly on the external field, if the latter is weak. In the most general 
case, we have 


D(x, t) = fax f ateo —x’,t—?')E(x’, t), (6) 


where we have implicitly assumed that all length scales are significantly larger than the crystal 
lattice, ensuring homogeneity. Thence, the response function depends only on the differences 
between spatial and temporal coordinates [2, 8]. In Fourier space, the convolutions turn into 
multiplications and the fields are decomposed into individual plane-wave components of the 
wavevector q and angular frequency w. Thus, in the Fourier domain, Eq. (6) reads 


D(q, w) = e(q, w)E(q, w). (7) 


For notational convenience, we designate the Fourier-transformed quantities with the same 
symbol as the original while they differ in the dependent variables. The Fourier transform of 
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an arbitrary field F(r,t) is given by F(r,t) = f F(q,w)e4*-“dqdt where w, q represent the 
Fourier transform quantities. Hence, the Fourier transform of the divergence equations of D 
and E yields 


—iq ` D(q,w) = ATP ol) (8) 
—iq - E(q,w) = 477,,.(q,@). (9) 


In longitudinal oscillations, the electron displacement field is in the direction of q (Figure 1(b)), 
thus, q: D = gD and q- E = gE, where D(q,w) and E(q,w) refer to longitudinal fields. Com- 
bining Eqs. (7)-(9) yields 

Pext(4,®) 


Pild) = “e(q.@) ` (10) 


Interestingly enough, in the absence of external charges, p,,,(q,w) = 0, Eq. (10) states that non- 
zero amplitudes of charge oscillation exist, that is, p,,,(q,~) 4 0, under the condition 


€(q,w) = 0. (11) 


In other words, in the absence of any external perturbation, free collective charge oscillations exist 
with dispersion relation w(q) that satisfies condition (11). These are plasmon modes, and conse- 
quently, Eq. (11) is referred as plasmon condition. Furthermore, condition (11) leads to E = —4nP, 
revealing that at plasmon frequencies the electric field is a pure depolarization field [1, 2]. 


We note that due to their longitudinal nature, plasmon waves cannot couple to any transverse 
wave such as electromagnetic waves; as a result, volume plasmons cannot be excited by light. 
On the other hand, moving charged particles can be used for exciting plasmons. For instance, 
an electron beam passing through a thin metal excites plasmons by transferring part of its 
energy to the plasmon excitation. As a result, plasmons do not decay directly via electromag- 
netic radiation but only through energy transfer to electron-hole excitation (Landau damping) 
[2, 8, 14]. 


2.2. Dynamical dielectric function 


Based on the plasmon condition (11), the problem has been reduced in the calculation of the 
dynamical dielectric function ¢(q, w). Further investigation of ¢(q, w) reveals the plasmon 
dispersion relation as well as the Landau-damping regime, that is, where plasmons decay very 
fast exciting electron-hole pairs [8]. Classically, in the long wavelength limit, the dielectric 
response €(0, w) can be calculated in the context of the plasma model [1, 11]. Let us consider 
the plasma model of Eq. (4) subjected to a weak and harmonic time-varying external field 
D(t) = D(w)e"'; Eq. (4) is modified to read 


mx (t) + 47te?nx(t) = —eD(t). (12) 


Assuming also a harmonic in time electron displacement, that is, x(t) = x(w)e, the Fourier 
transform of Eq. (12) yields 
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(—ma* + 4re*n)x(q,w) = —eD(q,w). (13) 


Introducing Eq. (3) in Eq. (13) and using the relation (7), we derive the spatially local dielectric 
response 


wp (0)* 


€(0,w) =1—- 2 


; (14) 
where the plasma frequency wp(0) is defined in Eq. (5). Eq. (14) verifies that the plasmon 


condition (11) is satisfied at the plasma frequency. The dielectric function (14) coincides with 
the Drude model permittivity. 


Further investigation of the dynamical dielectric function can be performed using quantum 
mechanics. An explicit form of ¢(q,w) including screening effect has been evaluated in the 
context of the random phase approximation (RPA) [8, 12-14] and is given by 


E(q,@) = 1 — ve(q)xo (q0) (15) 


where v, (q) is the Fourier transform of the Coulomb potential and x9(q,w) is the polarizability 
function, known as Lindhard formula [8, 12-14]. The Coulomb potential in two and three 
dimensions, respectively, reads 


ne (2D) 
Iqléo 
Ve(q) = 2 (16) 
Atte 
lql £r ee 


where €, represents the background lattice dielectric constant of the system. 


In RPA approach, the dynamical conductivity o(q,w) reads [8] 
iwe? 
o= “ge ola): (17) 


revealing the fundamental relation between ¢(q,w) and o(q,w) that also depends on system 
dimensions; we have finally 


e(q.0) = 1 +i olqw). (18) 


In the random phase approximation, the most important effect of interactions is that they 
produce electronic screening, while the electron-electron interaction is neglected. The polariz- 
ability of a non-interacting electron gas is represented by Lindhard formula as follows: 


f (€k+q) — f (€x) 
qa = = hw — a — &) + ihn _ 
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where factor 2 is derived by spin degeneracy (summation over the two possible values of spin 
s = Î,}) [8, 13, 14]. The summation is over all the wavevectors k, V is the volume, ifn repre- 
sents a small imaginary number to be brought to zero after the summation, and ex is the kinetic 
energy for the wave vector k. The carrier distribution f is given by Fermi-Dirac distribution 


-1 
f(a) = (exp Blax- u) +1) , where u is the chemical potential and ß =1/kgT with 


Boltzmann’s constant denoted by kg and T is the absolute temperature. Equation (19) describes 
processes in which a particle in state k, which is occupied with probability f (ex), is scattered 
into state k + q, which is empty with probability 1 — f (x44). Eqs. (15)-(19) consist of the basic 
equations for a detailed investigation of charge density fluctuations and the screening effect, 
electron-hole pair excitation, and plasmons. With respect to condition (11), the roots of Eq. (15) 
determine the plasmon modes. Moreover, the poles of 7) account for electron-hole pair excita- 
tion defining the Plasmon-damping regime [12-14]. 


For an analytical investigation, we split the summation of Eq. (19) in two parts. We make an 
elementary change of variables k + q — —k, in the term that includes f (€,;g), and assume that 
the kinetic energy is symmetric with respect to the wavevector, that is, ék = €_,. Therefore, 
formula (19) yields 


2 flex) F(&) 
w) = 20 
Xo(a) V (= hz — (€k+q — &) >, ħz + (€k+q — €k) (20) 
where z = w + in. At zero temperature, the chemical potential is equal to Fermi energy, that is, 
u = Er [8, 11, 14], and the Fermi-Dirac distribution is reduced to Heaviside step function, thus, 
f(ek)lr-0 = O(Er — ex). The kinetic energy of each electron of mass m in state k is given by 


A? |k|? 
=. (21) 
hence 
-Ë ae 2k 22 
a — = 5— (la)? + 2k: q). (22) 


At zero temperature, because of the Heaviside step function, the only terms that survive in 
summation (20) are those with |k| < ke, where ke is the Fermi wavenumber and related to 


Fermi energy by equation (21) as kp = (2mE¢/h?)'/?. Subsequently, we obtain for the Lindhard 
formula 


4 Ek+q — & 


7 7 
V i (A) — (ekg — &) 


X0(q.0) = (23) 


Summation turns into integration by using v>, ki (...) — (20)? f d°k(...), hence 
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o 4 3 Ek+q — Ek 
rola) = Er i 


where the imaginary part in z guarantees the convergence of the integrals around the poles 
hw = +(€k+q — €k). The poles of 79 determine the Landau-damping regime where plasmons 


decay into electron-hole pairs excitation. In particular, the damping regime is a continuum 
bounded by the limit values of (€k+q — €k); k takes its maximum absolute value |k| = kp and 


the inner product takes the extreme values kr k. q = +kelq|. 


h h 
— ZL (q= 2kr) < w < = (9 4+ 2ke), (25) 


where q = |q|. The Landau-damping continuum (electron-hole excitation regime) is demon- 
strated in Figure 2 by the shaded area. 


Introducing relation (22) into Eq. (24) and changing to spherical coordinates (r,0,0), where 
r = |k| and @ are the angle between k and q, we obtain 


Qn Et 2x COS o) sin 0 
Xo(q,w) = T af dx x A E 7 (26) 
Ta mz (=1 (i a + X cos 0) 


where x = r/kp is a dimensionless variable and vp = hk¢/m is the Fermi velocity. In the non- 
static (w>>vrq) and long wavelength (q«kr) limits, we can expand the integral in a power 
series of q. Keeping up to q? orders, we evaluate integral (26) and set the imaginary part of z 
zero, that is, z = w. That leads to a third-order approximation polarizability function 


0 2kFf q 


Figure 2. Dispersion relation of plasmons in the bulk of three-dimensional solid (blue solid line) and in two-dimensional 
electron gas (dashed red curve) plasmons. The shaded region demonstrates the Landau-damping regime where plasmons 
decay to electron-hole pairs excitation. 
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ee 2 372 2 
ole) = PE (14S), (27) 


which, in turn, yields the dielectric function by using formula (14) and the three-dimensional 
Coulomb interaction (16), hence 


e(q > 0,@) =1 cop(0) (1 2 (=), (28) 


w2 5 


where vacuum is assumed as background (£4 = 1) and we use the relation 0 [1, 15] where n is 
the electron density. The result (28) is reduced to simple Drude dielectric function (14) for 
q=0. 


The plasmon condition (11) determines the q-dependent plasmon dispersion relation wp(q). 
Demanding ¢(q,w) = 0, Eq. (28) yields approximately 


2 
woro (1435 (25) ). (29) 


Interestingly enough, the leading term of plasma frequency (29) does not include any quantum 
quantity, such as vr, which appears as non-local correction in sub-leading terms. That reveals 
that plasmons in 3D metals are purely classical modes. Moreover, a gap, that is, w,(0), appears 


in the plasmon spectrum of three-dimensional metals. The plasmon dispersion relation (29) is 
shown in Figure 2. 


In the random phase approximation, the electrons do not scatter, that is, collision between 
electrons and crystal impurities is not taken into account. As a consequence, the dielectric 
function is calculated to be purely real; this is nevertheless an unphysical result as can be seen 
clearly at zero frequency where the dielectric function is not well defined, that is, ¢(q,0) = %. 
The problem is cured by introducing a relaxation time t in the denominator of the dielectric 
function as follows: 


w,(q) 


EMG) = 17 SOD 


(30) 


We can phenomenologically prove expression (30) by using the simple plasma model. In 
particular, we modify the equation of motion (12) to a damped-driven harmonic oscillator by 
assuming that the motion of electron is damped via collisions occurring with a characteristic 
frequency y = 1/T [2]; this approach immediately leads to the dielectric response (30). Typi- 
cally values of relaxation time t are of the order 10~'* s, at room temperature. The relaxation 
time is determined experimentally. In the presence of t, the dielectric function (15) is well 
defined at w = 0, where the real part of permittivity has a peak with width t! known as 
Drude peak. Furthermore, it can be shown that equation (30) satisfies the Kramers-Kronig 
relations (sum rules) [1, 14, 15]. 
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2.3. Surface plasmon polariton 


A new guided collective oscillation mode called surface plasmon arises in the presence of a 
boundary. Surface plasmon is a surface electromagnetic wave that propagates along an inter- 
face between a conductor (metal) and an insulator (dielectric). This guided mode couples to 
electromagnetic waves resulting in a polariton. Surface plasmon polaritons (SPPs) occur at 
frequencies close to but smaller than plasma frequency. These surface modes show exceptional 
properties for applications of nanophotonics, specifically they constitute a class of nano- 
photonics themselves, namely nanoplasmonics. The basic property is the subwavelength 
nature, that is, the wavelength of SPPs is smaller than electromagnetic radiation at the same 
frequency and in the same medium [2, 3, 9]. 


Let us consider a waveguide formed by a planar interface at z = 0 consisting of two semi-infinite 
nonmagnetic media (permeability u = 1) with dielectric functions £1 and €? as Figure 3a 
denotes. The dielectric functions are assumed to be local in space (non-q—dependent) and non- 
local in time (w dependence), hence £12 = €12(w). Assuming harmonic in time dependence in 
the form u(r,t) = u(r)e', the Maxwell equations (in CGS units) in the absence of external 
charges and currents read 


V. (€jE;) =0 VxEj = ikoH; (31) 
V. (H;) =0 V xHj = —iejkoE; (32) 
where kp = w/c is the free space wavenumber and the index j denotes the media as j = 1 for 


z <0 and j=2 for z > 0. Combining Eqs. (31) and (32), the fields are decoupled into two 
separate Helmholtz equations [2, 4] as 


[V? + ke] ee =0 (33) 


where r = (x,y,z). For simplicity, let us assume surface electromagnetic waves propagating 
along one direction, chosen to be the x direction (Figure 3b), and show no spatial variations 
in the perpendicular in-plane direction, hence ðu = 0. Under this assumption, we are seeking 


Dielectric 


Figure 3. A planar interface is formed between a metal and a dielectric where surface plasmon polaritons (SPPs) 
propagate in (a) three- and (b) two-dimensional representation. (b) A schematic illustration of the SPP field. 
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electromagnetic waves of the form w(x) = tp (z)e%*, where tp; = (E;,H,)" and q will be the 
plasmon propagation constant. Substituting the aforementioned ansatz into Helmholtz equa- 
tion (33), we obtain the guided electromagnetic modes equation [2] 


s+ eH) (HO) <0 (4) 


Surface waves are waves that have been trapped at the interface (z = 0) and decay exponen- 
tially away from it (Ye) ~ e-SiFlfor k; > 0). Consequently, propagating wave solutions along 


z is not desired. In turn, we derive the surface wave condition 


Kj = VG — ke ER. (35) 


In order to determine the spatial field profiles and the SPP dispersion relation, we need to find 
explicit expressions for each field component of E and H. This can be achieved by solving the 
curl equations (31) and (32), which naturally lead to two self-consistent set of coupled 
governing equations. Each set corresponds to one of the fundamental polarizations, namely 
transverse magnetic (TM) (p-polarized waves) and transverse electric (TE) (s-polarized waves), 
hence 


Transverse magnetic (TM) Transverse electric (TE) 
7 > ade 
Ez = -74-Hy He =p 
of id 
io H; = 
Ex =-~——H j iv 
ú ko E j Oz y z ko Oz 
d Ey 2 — kpej)Ey =0 
52 Biv -= - ke) Hiy = 0 dz? ” (a = koej)Ey 


We focus on transverse magnetic (TM) polarization, in which the magnetic field H is parallel to 
the interface. Since the planar interface extends along (x,y) plane, the TM fields read 
E; = (E, 0,E;z) and H; = (0,Hj,,0). Solving the TM equations for surface waves, we obtain for 
each half plane 


ea (0 =1) 
ae (36) 
Ay = Ayei*es* 
ik, Ay ck 
Ez ee ng 1g, X Kz 7 
es ehe (37) 
q,4ı iq, x „kız 
E; — — ph 
e ehe (38) 
z>0 (G=2 
G=2) (39) 


H, = Adh ekz 
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ik2 A? igx ,—kp 
ey eg! & 40 
ko€2 i j ( ) 
q@rA2 igox ,—koz 
E, = — 2 dht ek 41 
Z ko€2 : i ( ) 


where k; is related to q; by Eq. (35). The boundary conditions imply that the parallel to interface 
components of electric (Ex) and magnetic (H,) fields must be continuous. Accordingly, we 
demand Eqs. (36) = (39) and Eqs. (37) = (40) at z = 0, hence we find the system of equations 


øn — eh Al 
ki ins k2 ah & ) =u (42) 


El €2 


which has a solution only if the determinant is zero. As an outcome, we obtain the so-called 
surface plasmon polariton condition 


i + ee 0. (43) 
E€ €2 

Condition (43) states that the interface must consist of materials with opposite signed permit- 
tivities, since surface wave condition requires the real part of both kı and kz to be non-negative 
numbers. For that reason, interface between metals and dielectrics may support surface 
plasmons, since metals show negative permittivity at frequencies smaller than plasma fre- 
quency [2]. Furthermore, boundary conditions demand the continuity of the normal to the 
interface electric displacement (D; = ¢;E;,) yielding the continuity of the plasmon propagation 
constant q4 = q, = q [4]. In turn, by combining Eq. (35) with Eq. (43) we obtain the dispersion 
relation for the surface plasmon polariton 


@ EJES 
_wv 44 
q(0) == foe (44) 


where €12 are, in general, complex functions of w. For a metal-dielectric interface, it is more 
convenient to use the notation £1 = €q and £2 = €,, for dielectric and metal permittivity, respec- 
tively. In long wavelengths, the SPP wavenumber is close to the light line in dielectric, viz. 
q=ko./éa, and the waves are extended over many wavelengths into the dielectrics [2, 4]; these 
waves are known as Sommerfeld-Zenneck waves and share similarities with free surface 
electromagnetic modes [2]. On the other hand, at the limit q — %, Eq. (44) asymptotically leads 
to the condition 


Egt Em = 0 (45) 


indicating the nonretarded surface plasmon limit [4]. In the vicinity of the nonretarded limit, 
Eq. (35) yields kj~q >> ko. Furthermore, in the nonretarded limit the phase velocity vph = w/q is 
tending to zero unveiling the electrostatic nature characterized by the surface plasmon [2, 3]. 
As a result, at the same frequency vpn is much smaller than the speed of light and, thus, the SPP 
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wavelength (Asp) is always smaller than the light wavelength (Apn), that is, Asp < Aph, reveal- 
ing the subwavelength nature of surface plasmon polaritons [2, 4]. In addition, due to the fact 
that SPP phase velocity is always smaller than the phase velocity of propagating electromag- 
netic waves, SPPs cannot radiate and, hence, they are well-defined surface propagating elec- 
tromagnetic waves. Demanding q — © in the dielectric function (30), we find the so-called 
surface plasmon frequency Wsp, which is the upper frequency limit that SPPs occur 


2 
a ~ @p 


2 
1+ €4 7 VIF €a 


Wsp = (46) 


indicating that SPPs always occur at frequencies smaller than bulk plasmons. 


If we follow the same procedure for transverse electric polarized fields, in which the electric 
field is parallel to interface and the only non-zero electromagnetic field components are E},Hx, 
and H., we will find the condition kı +k, = 0 [2]. This condition is satisfied only for 
ky = k = 0 unveiling that s-polarized surface modes do not exist. Consequently, surface 
plasmon polaritons are always TM electromagnetic waves. 


Due to metallic losses, SPPs decay exponentially along the interface restricting the propagation 
length. Mathematically speaking, losses are described by the small imaginary part in the 
complex dielectric function of metal €m = —Ew — ie”, where €/,,,€/,, > 0. Consequently, the SPPs 
propagation constant (44) becomes complex, that is, q = q' + ig", where the imaginary part 
accounts for losses of SPPs energy. In turn, the effective propagation length L, which shows the 
rate of change of the energy attenuation of SPPs [2, 3], is determined by the imaginary part 
Im|q] as L~' = 2Im|q]. 


Gain materials rather than passive regular dielectrics have been used to reduce the losses in 
SPP propagation. Gain materials are characterized by a complex permittivity function, that is, 
€q = E€} = +ie4, with €',e, > 0, where £4 is a small number compared to £4 and accounts for 
gain. As a result, gain dielectric gives energy to the system counterbalancing the metal losses. 
We investigate the SPP dispersion relation (44) in the presence of gain and loss materials, and 
find an explicit formula for gain e4 where the SPP wavenumber is reduced to real function, 
resulting in lossless SPPs propagation. In addition, we find an upper limit that values of gain 
are allowed. In this critical gain, the purely real SPP propagation constant becomes purely 
imaginary, destroying the SPP modes. 


The dispersion relation (44) can also be written as q = konsp [3], where nsp is the plasmon 


Nsp = 4| ; 47 
o Ed + Em ( ) 


We are seeking for a gain 4 such that the effective index nsp becomes real. Substituting the 


effective refractive index given by 


complex function describing the dielectric and metal into Eq. (47), the function nsp is written in 
the ordinary complex form as [3] 
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Nsp = 5 + isgn(y) a (48) 


where sgn(y) is the discontinuous signum function [3] and 


a SA = el leal? 


|Eq + ef, |? 


(49) 


ehje ft e feal? 
y= E SS (50) 
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with |z.| denoting the norm of the complex number z,. The poles in Eqs. (49) and (50) corre- 
spond to the nonretarded surface plasmon limit (45). 


Considering the plasmon effective index nsp in Eq. (48) in the (x, y) plane, we observe that 
lossless SPP propagation (Im|n,,] = Im[q] = 0) is warranted when the conditions y = 0 and 
x > 0 are simultaneously satisfied. Let us point out that for y = 0 and x < 0, although the 
imaginary part in Eq. (48) vanishes due to the signum function, its real part becomes imagi- 
nary, that is, nsp = iy |x|, which does not correspond to propagation waves. Solving Eq. (50) for 
y = 0 with respect to gain £4 and avoiding the nonretarded limit (45), that is, eg A —Em, we 
obtain two exact solutions [3] as follows: 


2 2 
+ Det el 
vg lem Pe |g [28am |) | 51 
a Bes, lal wy 


Due to the fact that £4 is real, we read from Eq. (51) that [3]. 


Jem)? > 24E% (52) 


Using inequality (52), we read for the solution £€4+ of Eq. (51) that e”4+ > €}. This is a contradic- 
tion since the €% is defined to be smaller than £4. Thus, €44 does not correspond to a physically 
relevant gain. 


Solving, on the other hand, Eq. (49) for x > 0, with respect to the dielectric gain €}, we 
determine a critical value £, distinguishing the regimes of lossless and prohibited SPP propa- 


gation [3], namely 
a lemt —1 (53) 
c — ca et ey, y 


hence, Eq. (53) sets an upper limit in values of gain. The appearance of critical gain can be 
understood as follows: In Eq. (51) the gain c4- becomes equal to critical gain €e when 
Ed + Em = 0 [3], where the last item is the nonretarded limit where q — ». Specifically, the 
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surface plasmon exists when the metal is characterized by the Drude dielectric function of 
Eq. (30), €4— = €c at w = wsp, corresponding to a maximum frequency [3]. 


In order to represent the above theoretical findings, we use the dielectric function of Eq. (30) to 
calculate the SPP dispersion relation for an interface consisting of silver with w,(0) = 13.67 
PHz and y = 0.1018 PHz, and silica glass with £} = 1.69 and for gain €} = £4- determined by 
Eq. (51). We represent in Figure 3a the SPP dispersion relation of Eq. (44) for lossless case 
(£4 = £4—), where the lossless gain is denoted by the inset image in Figure 3a. We indicate the 
real and imaginary of normalized SSP dispersion q/k, (kp = w,/c), with respect to the normal- 
ized frequency w/w). We observe, indeed, that for w < wsp the imaginary part of q vanishes, 
whereas for w > wsp the SPPs wavenumber is purely imaginary. Subsequently, in the vicinity 
of w = wsp a phase transition from lossless to prohibited SPPs propagation is expected [3]. 


We also solve numerically the full system of Maxwell equations (31) and (32) in a two-dimen- 
sional space for transverse magnetic polarization. The numerical experiments have been 
performed by virtue of the multi-physics commercial software COMSOL and the frequency w 
is confined in the range [0.3w,,0.75w,| with the integration step Aw = 0.01w,. In the same 


range, the lossless gain is calculated by Eq. (51), to be [8 -10~°, 8 - 107°]. For the excitation of 
SPPs on the metallic surface, we use the near-field technique [2, 3, 9, 10]. For this purpose, a 
circular electromagnetic source of radius R = 20 nm has been located 100 nm above the metal- 
lic surface acting as a point source, since the wavelength A of EM waves is much larger, that is, 
A >> R [2, 3]. In Figure 4b, we demonstrate, in a log-linear scale, the propagation length L, 
with respect to œw, subject in lossless gain eg— (blue line and open circles). For the sake of 
comparison, we plot L(w) in the absence of gain (green line and filled circles). The solid lines 
represent the theoretical predictions obtained by the definition of L, whereas the circles 
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Figure 4. (a) The surface plasmon polariton (SPP) dispersion relation q(w) in the presence of a gain material with gain 
corresponds to lossless SPP propagation. Re{q] and Im|q] are indicated by blue and red lines, respectively. The horizontal 
dashed black line denotes the SPP frequency (wsp = 0.61w,) where an interchanging between Re|q] and Im{q] appears. 
The dotted magenta line indicates the light line in the dielectric. (Inset) Demonstration of the gain leads to lossless SPP 
propagation. (b) Theoretical (solid lines) and numerical (circles) prediction of SPP propagation length L in the presence 
(blue) and in the absence (green) of gain dielectric showing a phase transition that happens at wsp (vertical dashed black 
line). Deviations between theoretical and numerical predications for w > wsp correspond to quasi-bound EM modes. The 
kp = Wp /c is used as normalized unit of wavenumbers and Wy as normalized unit for frequencies. 


18 Nanoplasmonics - Fundamentals and Applications 


indicate numerical results. For the numerical calculations, the characteristic propagation 
length has been estimated by the inverse of the slope of the Log(I), where I is the magnetic 
intensity along the interface [2-4]. The black vertical dashed line denotes the SPP resonance 
frequency wsp in which the phase transition appears. The graphs in Figure 4b indicate that in 
the presence of the lossless gain, SPPs may travel for very long, practically infinite, distances. 
Approaching the resonance frequency wsp L decreases rapidly leading to a steep phase transi- 
tion on the SPPs propagation. The deviations between theoretical and numerical results in 
Figure 4 for frequencies near and greater than wsp are attributed to the fact that in the regime 
Wsp < W < Wy, there are quasi-bound EM modes [2, 3], where EM waves are evanescent along 
the metal-dielectric interface and radiate perpendicular to it. Consequently, the observed EM 
field for w > wsp corresponds to radiating modes [3]. 


3. Two-dimensional plasmons 


In this section, we investigate plasmons in a two-dimensional electron gas (2DEG), where the 
electron sea is free to move only in two dimensions, tightly confined in the third. The reduced 
dimensions of electron confinement and Coulomb interaction cause crucial differences in 
plasmons excitation spectrum. For instance, plasmon spectrum in a 2DEG is gapless in contrast 
with three-dimensional case [13]. For the sake of completeness, we first discuss briefly 
plasmons in a regular 2DEG characterized by the usual parabolic dispersion relation (21) for a 
two-dimensional wavevector k lies in the plane of 2DEG. Thence, we focus on plasmons in a 
quite special two-dimensional material, viz. graphene. Graphene is a gapless two-dimensional 
semi-metal with linear dispersion relation. The linear energy spectrum offers great opportu- 
nity to describe graphene with chiral Dirac Hamiltonian for massless spin-1/2 fermions 
[7, 8, 10]. Furthermore, graphene can be doped with several methods, such as chemical doping 
[7], by applying an external voltage [10], or with lithium intercalation [16]. The doping shifts 
the Fermi level toward the conduction bands making graphene a great metal. The advantage to 
describe graphene electronic properties with massless carriers Dirac equation leads to excep- 
tional optical and electronic properties, like very high electric conductivity and ultra- 
subwavelength plasmons [6-8, 10]. 


3.1. Dynamical dielectric function of 2D metals 


In order to determine the plasmon spectrum of a two-dimensional electron gas, first of all we 
calculate the dielectric function in the context of random phase approximation (15) with v4 
being the two-dimensional Coulomb interaction of Eq. (16). In the Lindhard formula (23), V 
and k denote a two-dimensional volume and wave-vector, respectively. First, we investigate a 
2DEG described by the parabolic dispersion relation (21). The electrons are assumed to occupy 
a single band ignoring interband transitions, that is, transitions to higher bands. Thus, there is 
no orbital degeneracy (g,=1) resulting in the two-dimensional Fermi wavenumber 
kp = /2nn, where n is the carrier (electrons) density [13, 17]. Turning summation (23) into 
integral by the substitution Ve (...) = (2m)? { d@k(...), we obtain the Lindhard formula 


in integral form 
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_ 4 2 Ek+q — Ek 
Xo(qw) = (2)? fa [k| (hz)? Ba (Ak = ey) (54) 


The singe particle excitation continuum is still defined by expression (25), since the kinetic 
energy is considered to have the same form as in 3D case, even though the 2D Fermi wave- 
number has been modified. Transforming to polar coordinate system (r,@) and using relation 
(22), integral (54) reads 


rs — (22) ) (se + xcos o) a 


fass f” a L + 2x cos O 


where x is a dimensionless variable defined as x = r/kr. Previously, since we are interested in 
long wavelength limit (q « kf), we expand the integrand of Eq. (55) around q = 0. Keeping up 
to first orders of q, integral (55) yields 


(56) 


where z — w by sending the imaginary part of z to zero. The dielectric function is determined 
by the formula of Eq. (15) for 2D Coulomb interaction of Eq. (16), hence 


27ne"q 
si 57 
elpo) =1- (57) 
The 2DEG plasmon dispersion relation is determined by Eq. (11) to be 
2rne?q 
2D = APAS Mh 58 
TORR (58) 


related with volume plasmons dispersion relation by wP (q) = wpy/q/2 In contrast to three- 
dimensional electron gas where plasmon spectrum is gapped, in two-dimensional case the 
plasmon frequency depends on ,/q making the plasmon spectrum gapless. In Figure 2, the 2D 
plasmon dispersion relation (58) is demonstrated together with three-dimensional case. Fur- 


thermore, it is worth pointing out the similarity between the plasmon dispersion relation of 
2DEG of Eq. (58) and SPP of Eq. (44), that is, both show ,/q dependence. 


Let us now investigate the most special two-dimensional electron gas, namely graphene. At 
the limit where the excitation energy is small compared to Er, the dispersion relation of 
graphene, viz. the relation between kinetic energy ej and momentum p = hk, is described by 
two linear bands as 


ex = shoe|k| (59) 


where s = +1 indicates the conduction (+1) and valence (-1) band, respectively, vp is the two- 
dimensional Fermi velocity which is constant for graphene and equal to vp= 10° m/s 
[7, 8, 10, 16, 18]. Because of valley degeneracy g, = 2, the Fermi momentum is modified to read 
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kp = \/2nn/g¢, = vnn [8, 18]. The Fermi energy, given by Er = hurkp, becomes zero in the 
absence of doping (n = 0). As a consequence, the Er crosses the point where the linear valence 
and conduction bands touch each other, namely at the Dirac point, giving rise to the semi- 
metal character of the undoped graphene [7, 15, 16, 18]. The Lindhard formula of Eq. (19) 
needs to be generalized to include both intra- and interband transitions (valley degeneracy) as 
well as the overlap of states, hence 


( w) = — 882 ~~ Fh) — f(a) F (k,k + ) (60) 
ANEW Saha hw = (Egg =) tity I 


where the factors g, = g, =2 account to spin and valley degeneracy, respectively. The 
1 1 

Lindhard formula has been modified to contain two extra summations (> X ) 
s=—1 s'=-1 


corresponding to valley degeneracy for the two bands of Eq. (59). In addition, the overlap of 
states function F.y(k,k+q) has been introduced and defined by Fy(k,k+q) = 
(1+ ss’ cos y)/2, where y is the angle between k and k + q vectors [5, 18]. The term cos w can 
be expressed in |k], |k + q| and 0 terms, and subsequently the overlap function is written as [8] 


1 , |k| + |q| cos 0 
Fs (k,k ! q) =53 (1 + SS | | a ) (61) 


In long wavelength limit, we approximately obtain 


lqjcos@  |q|’ sin2@ 
k+qļ|=[|k|| 1+ a : 62 
k+ ql x( T (62) 


In this limit, we obtain for the graphene dispersion relation (59) the general form 


s= 


Ik] + lal( cos 0+ 3 sin20) J. (63) 


Gag -q= shor ( jk] 


S 


In turn, the plasmon-damping regimes are determined by the poles of polarizability (60) by 
substituting expression (63). Due to the valley degeneracy, there are two damping regimes 
corresponding, respectively, to intraband (s = s’) 


w < vrq (64) 


and interband (s = —s’) 


vp(2kp — q) < w < Up(2ke + q). (65) 


electron-hole pair excitations [8] demonstrated in Figure 5 by shaded areas. 


Substituting the long wavelength limit expression (62) in the overlap function (61), the latter 
reads 
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0 2kr q 


Figure 5. Blue solid line indicates the dispersion relation of graphene plasmons (w;”). The shaded regimes represent the 


intra- and interband Landau damping where plasmon decays to electron-hole pairs excitation. 


2 


i= = sin*@ ~1 s= 73' (intraband) 
Ezg (k,k + q) =] 2 (66) 
a sin*@ 20 ss’ (interband) 


Equation (66) states that in long wavelength limit, the interband contribution can be neglected 
[5], hence, the Lindhard formula (60) is simplified to 


xo(q > 0,0) = hes a Ae a l. (67) 


| 
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As it has already been mentioned, in zero temperature limit, the Fermi-Dirac distribution f (q 
is simplified to Heaviside step function O(kr F |k|). In this limit, the second term in the right 
hand of Eq. (67) is always zero, since O(kp + |k|) = O(ke + |k + q|) = 1, which reflects that all 
states in the valence band are occupied. Making again the elementary transformation 
k + q > —k in the term of Eq. (67) that includes f (e¢, q), we obtain 


xola 0,0) =S rak (68) 
0 ’ = $ 
V kt (hz)? — ag ety 
Turning summation (68) into integral, we read 
= e 
Xo(q > 0,w) (69) 
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Transforming to polar coordinates for r = |k| and using relation (63), we obtain the integral 
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d0, (70) 
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where x =r/kr, q= |q| and ņn = 0 =z = w. In non-static (w>>vrq) and long wavelength 
(q<kr) limits, we expand the integrator of Eq. (69) in series of q. Keeping up to first power of 
q/kr, we obtain 


OF 2m 
Xo(q — 0,w ae fa xf (x cos 0 + oe sin20) dé. (71) 


The evaluation of integral (71) is trivial and leads to the polarizability function of graphene 


2 
x(q > 0,0) =—5-G. (72) 


(73) 


indicating that at low energies doped graphene is described by a Drude-type dielectric func- 
tion with plasma frequency depending straightforward on the doping amount, namely the 
Fermi energy level Er. The plasma frequency of graphene monolayer is determined by condi- 
tion (11) and reads 


Gr = 2e2Er 
— Py 


(74) 


indicating the q'/? dependence likewise plasmons at a regular 2DEG. The most important 
result is the presence of A in the denominator of Eq. (74), which reveals that plasmons in 
graphene are purely quantum modes, that is, there are no classical plasmons in doped 


1/4 


graphene. In addition, graphene plasmon frequency is proportional to n’/*, which is different 


from classical 2D plasmon behavior where w3? ~ 


n!/2 [7, 18]. This is a direct consequence of 
the quantum relativistic nature of graphene, since Fermi energy is defined differently in any 
case, namely Er ~ kp ~ n!/2 in graphene, whereas, Er ~ ke ~ n in 2DEG case. In Figure 3(a), 


we represent the plasmon dispersion relation in doped graphene. 


3.2. Graphene plasmonic metamaterial 


Multilayers of plasmonic materials have been used for designing metamaterials providing elec- 
tromagnetic propagation behavior not found under normal circumstances like negative refrac- 
tion and epsilon-near-zero (ENZ) [9, 19, 20]. The bottleneck in creating plasmonic devices with 
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any desirable characteristic has been the limitations of typical 3D solids in producing perfect 
interfaces for the confinement of electrons and the features of dielectric host. This may no longer 
be a critical issue. The advent of truly two-dimensional materials like graphene (a metal), 
transition-metal dichalcogenides (TMDC’s, semiconductors), and hexagonal boron nitride (hBN, 
an insulator) makes it possible to produce structures with atomic-level control of features in the 
direction perpendicular to the stacked layers [9, 21]. This is ushering a new era in manipulating 
the properties of plasmons and designing devices with extraordinary behavior. 


Here, we propose a systematic method for constructing epsilon-near-zero (ENZ) metamaterials 
by appropriate combination on 2D materials. The aforementioned metamaterials exhibit inter- 
esting properties like diffractionless EM wave propagation with no phase delay [9]. We show 
analytically that EM wave propagation through layered heterostructures can be tuned dynam- 
ically by controlling the operating frequency and the doping level of the 2D metallic layers. 
Specifically, we find that multilayers of a plasmonic 2D material embedded in a dielectric host 
exhibit a plasmonic Dirac point (PDP), namely a point in wavenumber space where two linear 
coexisting dispersion curves cross each other, which, in turn, leads to an effective ENZ behav- 
ior [9]. To prove the feasibility of this design, we investigate numerically EM wave propagation 
in periodic plasmonic structures consisting of 2D metallic layers lying on yz plane in the form 
of graphene, arranged periodically along the x axis and possessing surface conductivity o;. The 
layers are embedded in a uniaxial dielectric host in the form of TMDC or hBN multilayers of 
thickness d and with uniaxial relative permittivity tensor Z4 with diagonal components 
Ex É Ey = Ez. We explore the resulting linear, elliptical, and hyperbolic EM dispersion relations 
which produce ENZ effect, ordinary and negative diffraction, respectively. 


We solve the analytical problem under TM polarization, with the magnetic field parallel to the 
y direction which implies that there is no interaction of the electric field with €. We consider a 
magnetically inert (relative permeability u = 1) lossless host (€,,¢, E R). For monochromatic 
harmonic waves in time, the Maxwell equations lead to three equations connecting the com- 
ponents of the E and H fields. For the longitudinal component [9, 19], Ez- = (ing /ko€z)(0H,/0x) 
where 1) = y Ho/€£o is the free space impedance. Defining the vector of the transversal field 


components as tp = (Ex, Hy)" gives [9] 


0 anai 
am) Ke Ox Ez Òx 
0 ~X SZ x 

ix b= kong ex > p (75) 
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Assuming EM waves propagating along the z axis, viz. w(x,z) = tp(x)ei#, Eq. (75) leads to an 
eigenvalue problem for the wavenumber k, of the plasmons along z [9, 19]. The metallic 2D 
planes are assumed to carry a surface current J, = 0,E,, which acts as a boundary condition in 
the eigenvalue problem. Furthermore, infinite number of 2D metals are considered to be 
arranged periodically, along x axis, with structural period d. The magnetic field reads 
H, (x)e* for —d < x < 0 and HY (x)e** for 0 < x < d on either side of the metallic plane at 


x = 0, with boundary conditions H; (0) — H, (0) = osE.(0) and òH; (0) = òH; (0). Due to the 
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periodicity, we use Bloch theorem along x as ch (x) = H, (x — d)e*4, with Bloch wavenumber 


kx. As a result, we arrive at the dispersion relation [9, 19, 20]: 


F(ky,k,) = cos (kxd) — cosh(«d) + = sinh(xa) =) (76) 


where x? = (¢./éx)(k — kjéx) expresses the anisotropy of the host medium and € = —(iosm)/ 
koéz) coincides with the so-called “plasmonic thickness” which determines the SPP decay 
length [9, 19, 20]. In particular, € is twice the SPP penetration length and defines the maximum 
distance between two metallic layers where the plasmons are strongly interacting [9, 19, 20]. 
We point out that for lossless 2D metallic planes gs is purely imaginary and € is purely real (for 
€z E R). At the center of the first Brillouin zone ky = 0, the equation has a trivial solution [19] for 
K = 0 > k; = kov\/éx which corresponds to the propagation of x-polarized fields travelling into 
the host medium with refractive index \/é, without interacting with the 2D planes which are 
positioned along z axis [22]. Near the Brillouin zone center (k,/kg <1 and «~0) and under the 
assumption of a very dense grid (d — 0), we obtain kyd «1 and xd «1, we Taylor expand the 
dispersion equation (76) to second order in d, hence 


ke d 
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The approximate relation (77) is identical to that of an equivalent homogenized medium 
described by dispersion: k2/ee + K/ee = k [9, 21]. Subsequently, from a metamaterial point 
of view, the entire system is treated as a homogeneous anisotropic medium with effective 
relative permittivities given by 


d= 
E e į 100s =e E (78) 


eff o- 
Ex = emar 
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We read from Eq. (78) the capability to control the behavior of the overall structure along the z 
direction. For instance, the choice d = €,/(€, — €x) leads to an isotropic effective medium with 
ceff L eff [9] 

Z x 3 


For the lossless case (Im|&] = 0), we identify two interesting regimes, viz. the strong plasmon 
coupling for d < € and the weak plasmon coupling for d > é. In both cases, plasmons develop 
along z direction at the interfaces between the conducting planes and the dielectric host. In the 
strong coupling case (d < €), plasmons of adjacent interfaces interact strongly with each other. 
As a consequence, the shape of the supported band of Eq. (77), in the (k,,k,) plane, is hyper- 
bolic (dashed red line in Figure 6(a)) and the system behaves as a hyperbolic metamaterial 
[9, 19, 22] with eS >0, ect < 0. On the other hand, in the weak plasmon coupling (d > €), the 
interaction between plasmons of adjacent planes is very weak. In this case, the shape of the 
dispersion relation (77) on the (k,,k,) plane is an ellipse (dotted black line in Figure 6(a)) and 
the systems act as an ordinary anisotropic media with ee ec > 0 [9]. We note that in the case 
€ < 0 the system does not support plasmons and the supported bands are always ellipses [9]. 
When either the 2D medium (Re[o;] 4 0) or the host material is lossy, a similar separation 
holds by replacing € by Re[é). 
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The most interesting case is the linear dispersion, where k, is linearly dependent on k, and 
dk,./dk, is constant for a wide range of k; [9, 19]. When this condition holds, the spatial 
harmonics travel with the same group velocity into the effective medium [9, 19]. To engineer 
our structure to exhibit a close-to-linear dispersion relation, we inspect the approximate ver- 
sion of Eq. (77): a huge coefficient for ky will make kj on the right-hand-side insignificant; if 


č = d, the term proportional to K increases without bound yielding a linear relation between k, 
and k,. With this choice, os = —i(kodéz/n)), and substituting in the exact dispersion relation 
Eq. (76), we find that (k,,k,) = (0,ko,/éx) becomes a saddle point for the transcendental func- 
tion F(k,,k,) giving rise to the conditions for the appearance of two permitted bands, namely 
two lines on the (k,,k,) plane across which F(k,,k,) = 0. This argument connects a mathemat- 
ical feature, the saddle point of the dispersion relation, with a physical feature, the crossing 
point of the two coexisting linear dispersion curves, the plasmonic Dirac point [9] (solid blue 
line in Figure 6(a)). From a macroscopic point of view, the choice € = d makes the effective 
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Figure 6. (a) The three supported dispersion plasmonic bands in (k,,k,) plane: hyberbolic (dashed red), elliptical (dotted 
black), and linear (solid blue) where plasmonic Dirac point (PDP) appears. (b) Combinations of graphene doping u, and 
free-space operational wavelengths A leading to epsilon-near-zero (ENZ) behavior (PDP in dispersion relation) for several 
lattice periods d (in nm). (c) Real and (d) imaginary parts of the effective permittivity £f for the choice d = 20 nm (dashed 
line in (b)); dashed curves indicate the ENZ regime. 
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permittivity along the z direction vanish, as is evident from Eq. (78). As a result, the existence 
of a PDP makes the effective medium behave like an ENZ material in one direction (ec = 0). 


The plasmonic length € is, typically, restricted in few nanometers (€ < 100 nm). Regular 
dielectrics always present imperfections in nanoscales, hence, the use of regular materials as 
dielectric hosts is impractical. Furthermore, graphene usually exfoliates or grows up on other 
2D materials. Because of the aforementioned reasons, it is strongly recommended that the 
dielectric host to be also a 2D material with atomic scale control of the thickness d and no 
roughness. For instance, one could build a dielectric host by stacking 2D layers of materials 
molybdenum disulfide (MoS2) [23] with essentially perfect planarity, complementing the pla- 
narity of graphene. 


Substituting the graphene dielectric function (73) into formula (18), we calculate the two- 
dimensional Drude-type conductivity of graphene [6, 19, 21] 


___ ie’, 
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where u, is the tunable chemical potential equal to Fermi energy Er and 7 is the transport- 
scattering time of the electrons [6, 19] introduced in the same manner as in Eq. (30). In what 
follows, we use bulk MoS, which at THz frequencies is assumed lossless with a diagonal 
permittivity tensor of elements, €x = 3.5 (out of plane) and €, = €; = 13 (in plane) [23]. 


The optical losses of graphene are taken into account using t = 0.5 ps [19]. Since the optical 
properties of the under-investigated system can be controlled by tuning the doping amount, 
the operating frequency or the structural period, in Figure 6(b), we show proper combinations 
of u, and operational wavelength in free space A which lead to a PDP for several values of 
lattice density distances (d = Re[é] in nm) [9]. To illustrate, for a reasonable distance between 
successive graphene planes of d = 20 nm, the real (Figure 6(c)) and imaginary (Figure 6(d)) 
effective permittivity values that can be emulated by this specific graphene-MoS, architecture 
determine the device characteristics at different frequencies and graphene-doping levels. Pos- 
itive values of Reļ[e%] are relatively moderate and occur for larger frequencies and lower 
doping levels of graphene; on the other hand, Im{e®*] is relatively small in the ENZ region as 
indicated by a dashed line in both graphs [9]. On the other hand, losses become larger as 
Re[e*"] gets more negative. 


To examine the actual electromagnetic field distribution in our graphene-MoS, configuration, 
we simulate the EM wave propagation through two finite structures consisting of 40 and 100 
graphene planes with Re[é] = 20.8nm and for operational wavelength in vacuum A = 12m 
(f = 25 THz = 0.1 eV). In order to have a complete picture of the propagation properties, we 
excite the under-investigating structures with a 2D dipole magnetic source as well as with a 
TM plane wave source. In particular, the 40-layered structure is excited by a 2D magnetic 
dipole source, which is positioned close to one of its two interfaces and oriented parallel to 
them, denoted by a white dot in Figure 7(a)—(c). On the other hand, the 100-layered configura- 
tion is excited by a plane source, which is located below the multilayer and is rotated by 20° 
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Figure 7. Spatial distribution of the magnetic field (color map) of graphene-MoS, multilayer structure located between 
the blue dashed lines and embedded in MoS, background. In (a)-(c), the metamaterial consists of 40 graphene sheets and 
excited by a magnetic dipole (white dot). In (d)-(f), the structure is composed by 100 graphene layers and excited by a TM 
plane wave source located at y = 0 and rotated 20° with respect to the interface. (a), (d) d = Re[&] (ENZ behavior). (b), (e) 
d = 0.7Re[é] hyperbolic metamaterial. (c), (f) d = 1.5Re[é] elliptical medium, where Re[é] = 20.8 nm. Due to high reflec- 
tions in (d), (e), we observe pattern formation of stationary waves below the metamaterial. 


with respect to the interface; the blue arrow in Figure 7(d) indicates the direction of the 
incident wave. The normalized to one spatial distribution of the magnetic field value is shown 
in Figure 7 in color representation, where the volume containing the graphene multilayers is 
between the dashed blue lines. To minimize the reflections, the background region is filled 
with a medium of the same dielectric properties as MoS). In Figure 7(a and d), the system is in 
the critical case (d = Re[é]), where the waves propagate through the graphene sheets without 
dispersion as in an ENZ medium. In Figure 7(b and e), the interlayer distance is d = 0.7Re[é] 
(strong plasmon-coupling regime) and the system shows negative (anomalous) diffraction. In 
Figure 7(c and f) d =1.5Re[E] (weak plasmon-coupling regime) and the EM wave show 
ordinary diffraction through the graphene planes [9]. 


4. Conclusion 


In summary, we have studied volume and surface plasmons beyond the classical plasma 
model. In particular, we have described electronic excitations in solids, such as plasmons and 
their damping mechanism, viz. electron-hole pairs excitation, in the context of the quantum 
approach random phase approximation (RPA), a powerful self-consistent theory for determin- 
ing the dielectric function of solids including screening non-local effect. The dielectric function 
and, in turn, the plasmon dispersion relation have been calculated for a bulk metal, a two- 
dimensional electron gas (2DEG) and for graphene, the famous two-dimensional semi-metal. 
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The completely different dispersion relation between plasmon in three- and two-dimensional 
metals has been pointed out. Furthermore, we have highlighted the fundamental difference 
between plasmons in a regular 2DEG and in doped graphene, indicating that plasmons in 
graphene are purely quantum modes, in contrast to plasmons in 2DEG, which originate from 
classical laws. Moreover, the propagation properties of surface plasmon polariton (SPP), a 
guided collective oscillation mode, have been also investigated. For the completeness of our 
theoretical investigation, we have outlined two applications. First, we have examined SPPs 
properties along an interface consisting of a bulk metal and an active (gain) dielectric. We have 
found that there is a gain value for which the metallic losses have been completely eliminated 
resulting in lossless SPP propagation. Second, we have investigated a plasmonic metamaterial 
composed of doped graphene monolayers. We have shown that depending on operating 
frequency, doping amount, and interlayer distance between adjacent graphene layers, the 
wave propagation properties present epsilon-near-zero behavior, normal, and negative refrac- 
tion, providing a metamaterial with tunable optical properties. 
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